% This m-script performs the monte carlo for the case with unobserved 
% networks in Lewbel, Qu and Tang (2018)
% Date : April 29, 2019
clear; % clc;
S       = 200;
L       = 240; % # of groups in sample 
n       = 10; % group size(s)
pn      = 1; % [.25 .25 .25 .25]; % distr'n of network size
Ln      = mnrnd(L,pn,1); % draw network sizes, e.g., Ln(1) = # of groups with size n(1)
lambda  = .7; alpha = 1; beta = [1.5,2,0]'; gamma = [.9,0,.6]';% social effects
mu_e    =  0; sig_e = 1; % param of structural error distr'n
% d       =  2*length(beta)+2; % # of structural parameters
p       = .5; % link formation prob
singularInd = zeros(S,1); % rng(5, 'twister');
rng(5, 'v5uniform'); disp('*******************');
[n,L],
for s = 1:S, % s, tic;
    % simulate groups with different sizes n = 10:5:25
    % --- X : a 4-by-1 cell, each element in the cell is a n(i)*2*Ln(i) array
    % --- y : a 4-by-1 cell, each element in the cell is a n(i)*Ln(i) array
    singularInd(s) = 0; % count pathological cases
    for i = 1 : length(n),
        e    =  normrnd(mu_e, sig_e, n(i), Ln(i)); 
        el   =  1;
        while el < Ln(i) + 1,       
           X{i}(:,:,el) = [randsample([-1 1 2]',n(i),1), normrnd(0,1,n(i),1), normrnd(1,2,n(i),1)];       
           Gt           =  binornd(1,p,[n(i),n(i)]);  
           Gs           =  Gt-diag(diag(Gt),0); 
           % skip pathological networks with isolated individuals
           if min(min(sum(Gs,2)),min(sum(Gs,1)))==0, 
                singularInd(s) = singularInd(s) + 1; 
                % disp('Warning: isolated individuals!'); 
                continue; 
           end;
           G            =  Gs./repmat(sum(Gs,2),1,n(i));
           M            = (eye(n(i))-lambda*G)\eye(n(i));
           y{i}(:,el)   =  M*(alpha+X{i}(:,:,el)*beta+G*X{i}(:,:,el)*gamma+e(:,el));
           el           =  el + 1;
        end;
    end;

    % Step one: calculate the reduced-form regression estimates
    [Mu_k,Mu_0] = getReducedForm(y{i},X{i});   k = size(Mu_k,3);  n = size(Mu_k,1);

    % Step two: use the reduced-form estimates to construct equation for structural parameters 
    [a,b] = calculate_ab(Mu_k);   
    
    % % actual (a_k,b_k)
    % for ki = 1:k-1,
    %    act_ab(:,ki) = [beta(ki),beta(k);gamma(ki),gamma(k)]\[1;-lambda];
    % end;
    % a = act_ab(1,:); b = act_ab(2,:); a, b, % a = [.6667, .5000]; b = [-1.6667, -.6667];    
    % sum(a(1)*diag(Mu_k(:,:,1))+b(1)*diag(Mu_k(:,:,3))), 
    % sum(a(2)*diag(Mu_k(:,:,2))+b(2)*diag(Mu_k(:,:,3))),
    % ODG = ones(n,n)-eye(n); 
    % (a(1)*sum(sum(ODG.*Mu_k(:,:,1),2)) + b(1)*sum(sum(ODG.*Mu_k(:,:,3),2)))/(n^2-n),
    % (a(2)*sum(sum(ODG.*Mu_k(:,:,2),2)) + b(2)*sum(sum(ODG.*Mu_k(:,:,3),2)))/(n^2-n), 

    % step three: calculate the closed-form estimates
    mk_hat   = ones(1,n)*Mu_k(:,:,k)*ones(n,1)/n; 
    Phi_hat  = [0 a(1) 0 b(1) zeros(1,3); 1 zeros(1,3) a(1) 0 b(1);...
                zeros(1,2) a(2) b(2) zeros(1,3); 1 zeros(1,4) a(2) b(2); ...
                mk_hat 0 0 1 0 0 1; zeros(1,3) 1 zeros(1,3); zeros(1,5) 1 0];
    v_hat    = [1 0 1 0 mk_hat 0 0]';
    coef_est =  Phi_hat\v_hat;
    int_est  =  mean(Mu_0)*(1-coef_est(1));
    Est(:,s) = [coef_est; int_est]; % toc;       
    
end;
MSE = mean((Est - repmat([lambda;beta;gamma;alpha],1,S)).^2,2);
% empirical distr'n of estimtes from simulations
MIV     = mean(Est,2); % empirical avg of IV estimator
bias    = MIV - [lambda; beta; gamma;alpha];
StdDev  = sqrt(diag(cov(Est'))); % empirical s.d of IV estimator
disp('M.S.E.  Bias   Std.Dev');
[MSE, bias, StdDev],
figure(1); 
subplot(3,2,1); hist(Est(1,:),50); title(sprintf('\\lambda = %g',lambda));
subplot(3,2,2); hist(Est(2,:),50); title(sprintf('\\beta_1 = %g',beta(1)));
subplot(3,2,3); hist(Est(3,:),50); title(sprintf('\\beta_2 = %g',beta(2)));
subplot(3,2,4); hist(Est(5,:),50); title(sprintf('\\gamma_1 = %g',gamma(1)));
subplot(3,2,5); hist(Est(7,:),50); title(sprintf('\\gamma_3 = %g',gamma(3)));
subplot(3,2,6); hist(Est(8,:),50); title(sprintf('\\alpha = %g',alpha));
suplabel(sprintf('Group size: %g ; Sample size: %g',n,L)); 
saveas(1,sprintf('unobs_nwk_n%d_L%d.fig',n,L));
save(sprintf('mc_unobs_nwk_n_%g_L_%g.mat',n,L));


